BARRETO et al. 2022: Habitat loss and biodiversity gain (BIOCON-S-22-00746)
Model formulas
For each of our response variables we run three candidate models for each scale:
Model 0 <- ~1, absence of effect
Model 1 <- percentage of forest cover at 3 and 5km radius
Model1.2 <- non-linear forest cover at 3 and 5km radius (quadratic model)
Loading data and gathering into one datasets calculated in script ‘data_wrang_GIT.Rmd’, content: 1) “beetles”: raw data for species occurrence per landscape/site; 2) “hab.data”: data containing classes of habitat association; 3) “env.data”: environmental variables, forest cover percent at 3 and 5km radius buffer 4) “div_hab”: diversity (alpha, beta and gamma) between classification habitat association 5) “rc.df”: Beta RC calculations 6) “ab_hab”: mean and median abundance between groups of habitat association 7) “st.hab_ab”: abundance data between classification habitat association per site
| Landscape | gamma_FS | gamma_NFS | alpha_FS | alpha_NFS | betad_FS | betad_NFS | ab.land_FS | ab.land_NFS | brc.fsNA | brc.nfsNA | brc.abund.fsNA | brc.abund.nfsNA | fc_3km | fc_5km |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 148 | 11 | 9 | 3.75 | 2.25 | 7.25 | 6.75 | 66 | 40 | -0.2989286 | 0.1920000 | -0.0067210 | 0.4313647 | 38.52523 | 43.86337 |
| 215 | 19 | 7 | 5.62 | 2.00 | 13.38 | 5.00 | 174 | 51 | 0.0989286 | -0.3042857 | 0.2793534 | 0.4023071 | 30.15485 | 25.23160 |
| 263 | 16 | 13 | 6.00 | 4.50 | 10.00 | 8.50 | 393 | 214 | -0.5085714 | -0.0057143 | 0.6216216 | 0.4039754 | 10.10641 | 13.41180 |
| 266 | 11 | 8 | 3.25 | 1.38 | 7.75 | 6.62 | 70 | 12 | -0.2342857 | 0.1352381 | 0.1091091 | 0.1504362 | 30.87525 | 31.16090 |
| 282 | 18 | 10 | 6.62 | 4.00 | 11.38 | 6.00 | 220 | 211 | -0.2835714 | -0.1852381 | 0.2157157 | 0.6922637 | 15.26270 | 16.35308 |
| 291 | 9 | 3 | 3.71 | 1.71 | 5.29 | 1.29 | 150 | 49 | -0.3938095 | -0.5809524 | 0.0873731 | -0.4176081 | 45.74294 | 45.13213 |
| 317 | 12 | 4 | 5.00 | 1.38 | 7.00 | 2.62 | 152 | 56 | -0.4033333 | -0.5364286 | 0.1798584 | -0.4283212 | 48.75805 | 47.85138 |
| 329 | 21 | 13 | 8.75 | 5.50 | 12.25 | 7.50 | 373 | 258 | -0.0692593 | 0.1207407 | 0.9381882 | 0.4310382 | 14.73791 | 13.62588 |
| 333 | 15 | 14 | 4.62 | 4.50 | 10.38 | 9.50 | 159 | 218 | -0.1096429 | 0.1907143 | 0.6733877 | 0.6187259 | 13.71403 | 14.50960 |
| 335 | 15 | 8 | 5.75 | 2.38 | 9.25 | 5.62 | 140 | 113 | -0.2660714 | -0.3266667 | 0.3353353 | -0.9401783 | 23.87913 | 19.21670 |
| 359 | 17 | 16 | 6.25 | 5.25 | 10.75 | 10.75 | 152 | 130 | -0.1721429 | 0.3191667 | 0.7743815 | 0.8145288 | 18.10208 | 17.62431 |
| 399 | 21 | 13 | 8.00 | 5.00 | 13.00 | 8.00 | 433 | 168 | -0.0852381 | 0.2525926 | 0.9580056 | 0.7037752 | 29.73616 | 24.97851 |
RESULTS
We found 28 species to be forest specialists, those that come from the Atlantic Forest domain and exclusively collected in forested habitats; while 23 were non-forest specialists, species from a broader domain and/or less restrictive to forest habitats (hereafter FS and NFS, respectively).
FOREST SPECIALISTS
Gamma diversity
Modelling
g0 <- glm(gamma_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(g.fs.aic <- AICctab(g0, g5k, g5kq, g3k, g3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## g5k -28.2 65.4 5.1 0.0 3 0.751
## g3k -30.0 68.9 3.3 3.5 3 0.127
## g5kq -28.1 69.9 5.2 4.6 4 0.076
## g0 -33.3 71.9 0.0 6.6 2 0.028
## g3kq -29.6 72.9 3.7 7.6 4 0.017
Summary top model
summary(g5k)##
## Call:
## glm(formula = gamma_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.24788 -0.14644 -0.02457 0.14244 0.30719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.46708 2.02439 10.604 9.26e-07 ***
## fc_5km -0.23174 0.05885 -3.938 0.00278 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.03547038)
##
## Null deviance: 0.80342 on 11 degrees of freedom
## Residual deviance: 0.34475 on 10 degrees of freedom
## AIC: 62.361
##
## Number of Fisher Scoring iterations: 3
g.fs <- g5k# rename top model to save and plot belowResiduals
sim <- simulateResiduals(g5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1271, p-value = 0.672
## alternative hypothesis: two.sided
hnp::hnp(g5k)## Gamma model
boot::glm.diag.plots(g5k) ## for glm poisson or neg binomSpatial autocorrelation All Moran’s I tests with our response variables to verify spatial auto-correlation resulted in non-significant which means non spatially auto-correlated.
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = 0.013940, expected = -0.090909, sd = 0.057965, p-value =
## 0.07048
## alternative hypothesis: Distance-based autocorrelation
Prediction
Alpha diversity
Modelling
a0 <- glm(alpha_FS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(a.fs.aic <- AICctab(a0, a5k, a5kq,
a3k, a3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## a5k -19.5 48.0 3.0 0.0 3 0.538
## a3k -20.5 50.0 2.0 2.0 3 0.200
## a0 -22.4 50.2 0.0 2.2 2 0.175
## a5kq -19.2 52.1 3.2 4.1 4 0.068
## a3kq -20.5 54.7 2.0 6.7 4 0.019
Summary top model
summary(a5k)##
## Call:
## glm(formula = alpha_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.43968 -0.11810 -0.03523 0.09321 0.36169
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.56549 0.96276 7.858 1.38e-05 ***
## fc_5km -0.07518 0.02853 -2.635 0.025 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.06116683)
##
## Null deviance: 1.00974 on 11 degrees of freedom
## Residual deviance: 0.62036 on 10 degrees of freedom
## AIC: 44.989
##
## Number of Fisher Scoring iterations: 4
a.fs <- a5k# rename top model to save and plot belowResiduals
sim <- simulateResiduals(a5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1023, p-value = 0.672
## alternative hypothesis: two.sided
hnp::hnp(a5k)## Gamma model
boot::glm.diag.plots(a5k)Spatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = 0.0054053, expected = -0.0909091, sd = 0.0566974, p-value =
## 0.08937
## alternative hypothesis: Distance-based autocorrelation
Prediction
Beta diversity additive
Modelling
### LM Beta for regular model selection
b0 <- glm(betad_FS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_FS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_FS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_FS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_FS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(b.fs.aic <- AICctab(b0, b5k, b5kq,
b3k, b3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## b5k -22.5 54.0 5.6 0.0 3 0.710
## b5kq -21.8 57.3 6.3 3.3 4 0.135
## b3k -24.4 57.7 3.7 3.8 3 0.108
## b3kq -23.3 60.2 4.8 6.2 4 0.031
## b0 -28.1 61.5 0.0 7.6 2 0.016
Summary top model
summary(b5k)##
## Call:
## glm(formula = betad_FS ~ fc_5km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.24292 -0.15104 -0.01746 0.05238 0.31092
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.93904 1.28413 10.855 7.46e-07 ***
## fc_5km -0.15787 0.03686 -4.283 0.0016 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.03504058)
##
## Null deviance: 0.83647 on 11 degrees of freedom
## Residual deviance: 0.33078 on 10 degrees of freedom
## AIC: 50.972
##
## Number of Fisher Scoring iterations: 4
b.fs <- b5k # rename top model to save and plot belowResiduals
sim <- simulateResiduals(b5k)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1273, p-value = 0.624
## alternative hypothesis: two.sided
hnp::hnp(b5k) ## for glm Gamma or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(b5k) ## for glm Gamma or neg binomSpatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = 0.021701, expected = -0.090909, sd = 0.057512, p-value =
## 0.05023
## alternative hypothesis: Distance-based autocorrelation
Prediction
bRaup-Crick (presence/absence-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.fsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.fsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc0 4.8 -4.2 0.0 0.0 2 0.419
## rc3kq 8.4 -3.0 3.6 1.2 4 0.227
## rc5k 5.5 -2.0 0.7 2.2 3 0.139
## rc5kq 7.8 -1.9 3.0 2.4 4 0.129
## rc3k 5.0 -1.1 0.2 3.2 3 0.086
Summary of 2nd best model, as the top is of absence of effect (~1)
summary(rc3kq)##
## Call:
## lm(formula = brc.fsNA ~ fc_3km + I(fc_3km^2), data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.16680 -0.10052 -0.01060 0.07993 0.21297
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.6934193 0.2287612 -3.031 0.0142 *
## fc_3km 0.0426487 0.0178385 2.391 0.0405 *
## I(fc_3km^2) -0.0007772 0.0003009 -2.583 0.0296 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1391 on 9 degrees of freedom
## Multiple R-squared: 0.449, Adjusted R-squared: 0.3266
## F-statistic: 3.667 on 2 and 9 DF, p-value: 0.06842
rc1.fs <- rc0# rename top model to save and plot belowResiduals
sim <- simulateResiduals(rc3kq)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.82544, p-value = 0.816
## alternative hypothesis: two.sided
hnp::hnp(rc3kq)## Gaussian model (lm object)
Spatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = -0.149708, expected = -0.090909, sd = 0.057697, p-value =
## 0.3082
## alternative hypothesis: Distance-based autocorrelation
Prediction
Total abundance
Modelling
### LM "gamma" abundanec for regular model selection
gab0 <- glm.nb(ab.land_FS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_FS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_FS ~ fc_5km + I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_FS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_FS ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(gab.fs.aic <- AICctab(gab0, gab5k, gab5kq,
gab3k, gab3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## gab5k -70.6 150.1 2.1 0.0 3 0.404
## gab0 -72.7 150.7 0.0 0.5 2 0.311
## gab3k -71.2 151.4 1.5 1.2 3 0.218
## gab5kq -70.4 154.6 2.2 4.5 4 0.044
## gab3kq -71.1 155.8 1.6 5.7 4 0.023
Summary 2nd best model, as top model was the one of absence of effect (~1):
summary(gab5k)##
## Call:
## glm.nb(formula = ab.land_FS ~ fc_5km, data = data_div, init.theta = 4.550032026,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6556 -0.9899 -0.2919 0.6493 1.8325
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.93867 0.31886 18.62 <2e-16 ***
## fc_5km -0.02505 0.01108 -2.26 0.0238 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.55) family taken to be 1)
##
## Null deviance: 17.357 on 11 degrees of freedom
## Residual deviance: 12.410 on 10 degrees of freedom
## AIC: 147.13
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.55
## Std. Err.: 1.84
##
## 2 x log-likelihood: -141.134
gab.fs <- gab5k# rename top model to save and plot belowResiduals
sim <- simulateResiduals(gab5k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 1.1526, p-value = 0.496
## alternative hypothesis: two.sided
hnp::hnp(gab5k)## Negative binomial model (using MASS package)
boot::glm.diag.plots(gab5k)Spatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = -0.075229, expected = -0.090909, sd = 0.058432, p-value =
## 0.7884
## alternative hypothesis: Distance-based autocorrelation
# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab3k.r2 <- r2(gab3k)
## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance FS",
"FC at 5km",
round(gab.fs.aic$AICc[1],2),
round(gab.fs.aic$dAICc[1],1),
gab.fs.aic$df[1],
round(gab.fs.aic$weight[1],2),
round(gab5k.r2$R2_Nagelkerke,2)),
c(" ", "~ 1 (NULL)",
round(gab.fs.aic$AICc[2],2),
round(gab.fs.aic$dAICc[2],1),
gab.fs.aic$df[2],
round(gab.fs.aic$weight[2],2),
"-"),
c(" ", "FC at 3km",
round(gab.fs.aic$AICc[3],2),
round(gab.fs.aic$dAICc[3],1),
gab.fs.aic$df[3],
round(gab.fs.aic$weight[3],2),
round(gab3k.r2$R2_Nagelkerke,2)))Prediction
bRaup-Crick (abundance-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.abund.fsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.fsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.fsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.fsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.fsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.fs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc5k 0.1 8.8 3.8 0.0 3 0.595
## rc3k -0.8 10.7 2.9 1.8 3 0.237
## rc0 -3.7 12.7 0.0 3.9 2 0.083
## rc5kq 0.2 13.3 3.9 4.5 4 0.062
## rc3kq -0.8 15.4 2.9 6.6 4 0.022
Summary top model:
summary(rc5k)##
## Call:
## lm(formula = brc.abund.fsNA ~ fc_5km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39189 -0.17995 -0.01784 0.15689 0.50748
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.905410 0.176830 5.120 0.000451 ***
## fc_5km -0.018211 0.006124 -2.974 0.013955 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2629 on 10 degrees of freedom
## Multiple R-squared: 0.4693, Adjusted R-squared: 0.4162
## F-statistic: 8.843 on 1 and 10 DF, p-value: 0.01396
rc2.fs <- rc5k# rename top model to save and plot belowResiduals
sim <- simulateResiduals(rc5k)
plot(sim)testDispersion(sim) # ok##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc5k)## Gaussian model (lm object)
Spatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = -0.076074, expected = -0.090909, sd = 0.057165, p-value =
## 0.7952
## alternative hypothesis: Distance-based autocorrelation
Prediction
NON-FOREST SPECIALISTS
Gamma diversity
Modelling
### LM alpha for regular model selection
g0 <- glm(gamma_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
g5k <- glm(gamma_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
g5kq <- glm(gamma_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
g3k <- glm(gamma_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
g3kq <- glm(gamma_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(g.nfs.aic <- AICctab(g0,
g5k, g5kq,
g3k, g3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## g3k -26.3 61.6 7.6 0.0 3 0.7974
## g5k -28.4 65.9 5.5 4.2 3 0.0958
## g3kq -26.1 65.9 7.8 4.3 4 0.0950
## g5kq -28.4 70.5 5.5 8.9 4 0.0093
## g0 -33.9 73.2 0.0 11.6 2 0.0025
Summary top model:
summary(g3k)##
## Call:
## glm(formula = gamma_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.38295 -0.23805 -0.05072 0.13490 0.38893
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.39301 1.95020 8.919 4.49e-06 ***
## fc_3km -0.28151 0.04841 -5.815 0.000169 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.07298721)
##
## Null deviance: 2.46666 on 11 degrees of freedom
## Residual deviance: 0.71058 on 10 degrees of freedom
## AIC: 58.646
##
## Number of Fisher Scoring iterations: 5
g.nfs <- g3k # rename top model to save and plot belowResiduals
sim <- simulateResiduals(g3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.89592, p-value = 0.936
## alternative hypothesis: two.sided
hnp::hnp(g3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(g3k) ## for glm poisson or neg binomSpatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = -0.123318, expected = -0.090909, sd = 0.057839, p-value =
## 0.5753
## alternative hypothesis: Distance-based autocorrelation
Prediction
Alpha diversity
Modelling
### LM alpha for regular model selection
a0 <- glm(alpha_NFS ~ 1, data=data_div, family= Gamma (link = "identity"))
a5k <- glm(alpha_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
a5kq <- glm(alpha_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
a3k <- glm(alpha_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
a3kq <- glm(alpha_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(a.nfs.aic <- AICctab(a0, a5k, a5kq,
a3k, a3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## a3k -16.1 41.2 5.7 0.0 3 0.497
## a5k -16.4 41.8 5.4 0.6 3 0.365
## a5kq -15.7 45.1 6.1 3.9 4 0.070
## a3kq -15.9 45.5 5.9 4.3 4 0.057
## a0 -21.8 49.0 0.0 7.8 2 0.010
Summary top model:
summary(a3k)##
## Call:
## glm(formula = alpha_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.66790 -0.14831 0.00666 0.10705 0.54796
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.65928 0.79667 7.104 3.28e-05 ***
## fc_3km -0.08855 0.02018 -4.388 0.00136 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1090929)
##
## Null deviance: 2.9678 on 11 degrees of freedom
## Residual deviance: 1.1728 on 10 degrees of freedom
## AIC: 38.181
##
## Number of Fisher Scoring iterations: 5
a.nfs <- a3k # rename top model to save and plot belowResiduals
sim <- simulateResiduals(a3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.94044, p-value = 0.88
## alternative hypothesis: two.sided
hnp::hnp(a3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(a3k) ## for glm poisson or neg binomSpatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = -0.099144, expected = -0.090909, sd = 0.056241, p-value =
## 0.8836
## alternative hypothesis: Distance-based autocorrelation
Prediction
Beta diversity additive
Modelling
### LM Beta for regular model selection
b0 <- glm(betad_NFS~ 1, data=data_div, family= Gamma (link = "identity"))
b5k <- glm(betad_NFS ~ fc_5km, data=data_div, family= Gamma (link = "identity"))
b5kq <- glm(betad_NFS ~ fc_5km + I(fc_5km^2), data=data_div, family= Gamma (link = "identity"))
b3k <- glm(betad_NFS ~ fc_3km, data=data_div, family= Gamma (link = "identity"))
b3kq <- glm(betad_NFS ~ fc_3km + I(fc_3km^2), data=data_div, family= Gamma (link = "identity"))
## AIC-based model selection
(b.nfs.aic <- AICctab(b0, b5k, b5kq,
b3k, b3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## b3k -24.4 57.7 5.5 0.0 3 0.734
## b3kq -23.8 61.4 6.0 3.7 4 0.118
## b5k -26.2 61.4 3.7 3.7 3 0.115
## b0 -29.9 65.1 0.0 7.4 2 0.019
## b5kq -26.0 65.6 3.9 7.9 4 0.014
Summary top model:
summary(b3k)##
## Call:
## glm(formula = betad_NFS ~ fc_3km, family = Gamma(link = "identity"),
## data = data_div)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.71087 -0.18354 -0.04943 0.17316 0.48985
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 11.78046 1.59008 7.409 2.29e-05 ***
## fc_3km -0.19439 0.03898 -4.987 0.000548 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1092755)
##
## Null deviance: 2.9582 on 11 degrees of freedom
## Residual deviance: 1.2094 on 10 degrees of freedom
## AIC: 54.73
##
## Number of Fisher Scoring iterations: 6
b.nfs <- b3k # rename top model to save and plot belowResiduals
sim <- simulateResiduals(b3k)
plot(sim)testDispersion(sim)##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.65169, p-value = 0.624
## alternative hypothesis: two.sided
hnp::hnp(b3k) ## for glm poisson or neg binom (maybe also others)## Gamma model
boot::glm.diag.plots(b3k) ## for glm poisson or neg binomSpatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = -0.146700, expected = -0.090909, sd = 0.057505, p-value =
## 0.332
## alternative hypothesis: Distance-based autocorrelation
Prediction
bRaup-Crick (presence/absence-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.nfsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.nfsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc3k -0.6 10.1 2.0 0.0 3 0.361
## rc0 -2.6 10.6 0.0 0.4 2 0.293
## rc5k -1.1 11.3 1.5 1.1 3 0.207
## rc3kq 0.6 12.5 3.2 2.4 4 0.110
## rc5kq -0.8 15.3 1.8 5.1 4 0.028
Summary of 2nd best model, 1st wqs of absence of effect (~1):
summary(rc3k)##
## Call:
## lm(formula = brc.nfsNA ~ fc_3km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.3020 -0.2156 -0.0678 0.2557 0.4084
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.287950 0.190935 1.508 0.1625
## fc_3km -0.013092 0.006505 -2.013 0.0719 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2781 on 10 degrees of freedom
## Multiple R-squared: 0.2883, Adjusted R-squared: 0.2171
## F-statistic: 4.051 on 1 and 10 DF, p-value: 0.07185
rc1.nfs <- rc3k # rename top model to save and plot belowResiduals
sim <- simulateResiduals(rc3k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)## Gaussian model (lm object)
Spatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = -0.159951, expected = -0.090909, sd = 0.058625, p-value =
## 0.2389
## alternative hypothesis: Distance-based autocorrelation
Prediction
Total abundance
Modelling
### LM "gamma" abundanec for regular model selection
gab0 <- glm.nb(ab.land_NFS ~ 1, data=data_div)
gab5k <- glm.nb(ab.land_NFS ~ fc_5km, data=data_div)
gab5kq <- glm.nb(ab.land_NFS ~ fc_5km + I(fc_5km^2), data=data_div)
gab3k <- glm.nb(ab.land_NFS ~ fc_3km, data=data_div)
gab3kq <- glm.nb(ab.land_NFS ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(gab.nfs.aic <- AICctab(gab0, gab5k, gab5kq,
gab3k, gab3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## gab5k -63.9 136.7 5.3 0.0 3 0.402
## gab3k -64.2 137.4 4.9 0.7 3 0.285
## gab5kq -62.0 137.7 7.1 0.9 4 0.253
## gab3kq -63.6 141.0 5.5 4.3 4 0.047
## gab0 -69.1 143.6 0.0 6.8 2 0.013
Summary top model:
summary(gab5k)##
## Call:
## glm.nb(formula = ab.land_NFS ~ fc_5km, data = data_div, init.theta = 4.180831522,
## link = log)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.8524 -0.4383 0.2265 0.4858 0.8634
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.92551 0.33598 17.64 < 2e-16 ***
## fc_5km -0.04812 0.01180 -4.08 4.5e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(4.1808) family taken to be 1)
##
## Null deviance: 29.171 on 11 degrees of freedom
## Residual deviance: 12.765 on 10 degrees of freedom
## AIC: 133.73
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 4.18
## Std. Err.: 1.78
##
## 2 x log-likelihood: -127.727
gab.nfs <- gab5k # rename top model to save and plot belowResiduals
sim <- simulateResiduals(gab5k)
plot(sim)testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.43345, p-value = 0.192
## alternative hypothesis: two.sided
hnp::hnp(gab5k)## Negative binomial model (using MASS package)
boot::glm.diag.plots(gab5k)Spatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = -0.079003, expected = -0.090909, sd = 0.055904, p-value =
## 0.8313
## alternative hypothesis: Distance-based autocorrelation
# R2 based on Nakagawa et al. 2017: variação explicada pelos efeitos fixos dividida pela total
gab5k.r2 <- r2(gab5k)
gab5kq.r2 <- r2(gab5kq)
gab3k.r2 <- r2(gab3k)
## build AIC table do add model selection results:
aic.tab <- rbind(aic.tab, c("Total abundance NFS",
"FC at 5km",
round(gab.nfs.aic$AICc[1],2),
round(gab.nfs.aic$dAICc[1],1),
gab.nfs.aic$df[1],
round(gab.nfs.aic$weight[1],2),
round(gab5k.r2$R2_Nagelkerke,2)),
c(" ", "FC at 3km",
round(gab.nfs.aic$AICc[2],2),
round(gab.nfs.aic$dAICc[2],1),
gab.nfs.aic$df[2],
round(gab.nfs.aic$weight[2],2),
round(gab3k.r2$R2_Nagelkerke,2)),
c(" ", "Non-linear FC at 5km",
round(gab.nfs.aic$AICc[3],2),
round(gab.nfs.aic$dAICc[3],1),
gab.nfs.aic$df[3],
round(gab.nfs.aic$weight[3],2),
round(gab5kq.r2$R2_Nagelkerke,2)
))Prediction
bRaup-Crick (abundance-based)
Modelling
### LM gamma for regular model selection
rc0 <- lm(brc.abund.nfsNA ~ 1, data=data_div)
rc5k <- lm(brc.abund.nfsNA ~ fc_5km, data=data_div)
rc5kq <- lm(brc.abund.nfsNA ~ fc_5km + I(fc_5km^2), data=data_div)
rc3k <- lm(brc.abund.nfsNA ~ fc_3km, data=data_div)
rc3kq <- lm(brc.abund.nfsNA ~ fc_3km + I(fc_3km^2), data=data_div)
## AIC-based model selection
(rc.nfs.aic <- AICctab(rc0, rc5k, rc5kq,
rc3k, rc3kq,
base=T,weights=T, logLik=T))## logLik AICc dLogLik dAICc df weight
## rc3k -7.4 23.8 1.9 0.0 3 0.367
## rc0 -9.3 23.9 0.0 0.1 2 0.356
## rc5k -8.0 24.9 1.3 1.1 3 0.215
## rc3kq -7.3 28.3 2.0 4.5 4 0.040
## rc5kq -7.9 29.5 1.4 5.6 4 0.022
Summary 2nd best model, 1st is of absence of effect (~1):
summary(rc3k)##
## Call:
## lm(formula = brc.abund.nfsNA ~ fc_3km, data = data_div)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.23916 -0.18519 0.05082 0.27801 0.53338
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.82320 0.33770 2.438 0.0350 *
## fc_3km -0.02195 0.01151 -1.908 0.0855 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4918 on 10 degrees of freedom
## Multiple R-squared: 0.2669, Adjusted R-squared: 0.1936
## F-statistic: 3.641 on 1 and 10 DF, p-value: 0.08547
rc2.nfs <- rc3k # rename top model to save and plot belowResiduals
sim <- simulateResiduals(rc3k)
plot(sim)## qu = 0.75, log(sigma) = -3.631427 : outer Newton did not converge fully.
testDispersion(sim) # a lil under##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## dispersion = 0.91715, p-value = 0.984
## alternative hypothesis: two.sided
hnp::hnp(rc3k)## Gaussian model (lm object)
Spatial autocorrelation
testSpatialAutocorrelation(sim, x = env.data$x, y = env.data$y)##
## DHARMa Moran's I test for distance-based autocorrelation
##
## data: sim
## observed = -0.118475, expected = -0.090909, sd = 0.055483, p-value =
## 0.6193
## alternative hypothesis: Distance-based autocorrelation
Prediction
AICc-model selection table
| Diversity | Model | AICc | dAICc | df | weight | r2 |
|---|---|---|---|---|---|---|
| Gamma FS | FC at 5km | 65.36 | 0 | 3 | 0.75 | 0.58 |
| Alpha FS | FC at 5km | 47.99 | 0 | 3 | 0.54 | 0.4 |
| FC at 3km | 49.97 | 2 | 3 | 0.2 | 0.29 | |
| Beta FS | FC at 5km | 53.97 | 0 | 3 | 0.71 | 0.61 |
| RCbeta (P/A-based) FS | ~ 1 (NULL) | -4.25 | 0 | 2 | 0.42 |
|
| Non-linear FC at 3km | -3.02 | 1.2 | 4 | 0.23 | 0.33 | |
| Total abundance FS | FC at 5km | 150.13 | 0 | 3 | 0.4 | 0.44 |
| ~ 1 (NULL) | 150.66 | 0.5 | 2 | 0.31 |
|
|
| FC at 3km | 151.37 | 1.2 | 3 | 0.22 | 0.33 | |
| RCbeta (abund-based) FS | FC at 5km | 8.81 | 0 | 3 | 0.6 | 0.42 |
| FC at 3km | 10.65 | 1.8 | 3 | 0.24 | 0.32 | |
| Gamma NFS | FC at 3km | 61.65 | 0 | 3 | 0.8 | 0.73 |
| Alpha NFS | FC at 3km | 41.18 | 0 | 3 | 0.5 | 0.63 |
| FC at 5km | 41.8 | 0.6 | 3 | 0.37 | 0.61 | |
| Beta NFS | FC at 3km | 57.73 | 0 | 3 | 0.73 | 0.62 |
| RCbeta (P/A-based) NFS | FC at 3km | 10.15 | 0 | 3 | 0.36 | 0.22 |
| ~ 1 (NULL) | 10.56 | 0.4 | 2 | 0.29 |
|
|
| FC at 5km | 11.26 | 1.1 | 3 | 0.21 | 0.14 | |
| Total abundance NFS | FC at 5km | 136.73 | 0 | 3 | 0.4 | 0.82 |
| FC at 3km | 137.42 | 0.7 | 3 | 0.28 | 0.79 | |
| Non-linear FC at 5km | 137.66 | 0.9 | 4 | 0.25 | 0.93 | |
| RCbeta (abund-based) NFS | FC at 3km | 23.83 | 0 | 3 | 0.37 | 0.19 |
| ~ 1 (NULL) | 23.89 | 0.1 | 2 | 0.36 |
|
|
| FC at 5km | 24.9 | 1.07 | 3 | 0.22 | 0.12 |
MODEL COEFFICIENTS TABLE
# fazer tabela de coeficientes dos modelos selectionados
(tab_gfs <- tab_model(g.fs, show.p = F, show.se = TRUE))| gamma FS | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | 21.47 | 2.02 | 17.69 – 25.70 |
| fc 5km | -0.23 | 0.06 | -0.35 – -0.11 |
| Observations | 12 | ||
| R2 Nagelkerke | 0.579 | ||
(tab_afs <- tab_model(a.fs, show.p = F, show.se = TRUE))| alpha FS | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | 7.57 | 0.96 | 5.84 – 9.57 |
| fc 5km | -0.08 | 0.03 | -0.13 – -0.02 |
| Observations | 12 | ||
| R2 Nagelkerke | 0.396 | ||
(tab_bfs <- tab_model(b.fs, show.p = F, show.se = TRUE))| betad FS | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | 13.94 | 1.28 | 11.50 – 16.66 |
| fc 5km | -0.16 | 0.04 | -0.23 – -0.08 |
| Observations | 12 | ||
| R2 Nagelkerke | 0.613 | ||
(tab_brc1fs <- tab_model(rc1.fs, show.p = F, show.se = TRUE))| brc fs NA | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | -0.23 | 0.05 | -0.33 – -0.12 |
| Observations | 12 | ||
| R2 / R2 adjusted | 0.000 / 0.000 | ||
(tab_gabfs <- tab_model(gab.fs, show.p = F, show.se = TRUE))| ab land FS | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | std. Error | CI |
| (Intercept) | 379.43 | 120.98 | 208.54 – 714.06 |
| fc 5km | 0.98 | 0.01 | 0.96 – 1.00 |
| Observations | 12 | ||
| R2 Nagelkerke | 0.442 | ||
(tab_brc2fs <- tab_model(rc2.fs, show.p = F, show.se = TRUE))| brc abund fs NA | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | 0.91 | 0.18 | 0.51 – 1.30 |
| fc 5km | -0.02 | 0.01 | -0.03 – -0.00 |
| Observations | 12 | ||
| R2 / R2 adjusted | 0.469 / 0.416 | ||
# Nnfs
(tab_gnfs <- tab_model(g.nfs, show.p = F, show.se = TRUE))| gamma NFS | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | 17.39 | 1.95 | 13.79 – 21.65 |
| fc 3km | -0.28 | 0.05 | -0.38 – -0.18 |
| Observations | 12 | ||
| R2 Nagelkerke | 0.733 | ||
(tab_anfs <- tab_model(a.nfs, show.p = F, show.se = TRUE))| alpha NFS | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | 5.66 | 0.80 | 4.28 – 7.39 |
| fc 3km | -0.09 | 0.02 | -0.13 – -0.05 |
| Observations | 12 | ||
| R2 Nagelkerke | 0.634 | ||
(tab_bnfs <- tab_model(b.nfs, show.p = F, show.se = TRUE))| betad NFS | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | 11.78 | 1.59 | 8.81 – 15.36 |
| fc 3km | -0.19 | 0.04 | -0.28 – -0.11 |
| Observations | 12 | ||
| R2 Nagelkerke | 0.621 | ||
(tab_brc1nfs <- tab_model(rc1.nfs, show.p = F, show.se = TRUE))| brc nfs NA | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | 0.29 | 0.19 | -0.14 – 0.71 |
| fc 3km | -0.01 | 0.01 | -0.03 – 0.00 |
| Observations | 12 | ||
| R2 / R2 adjusted | 0.288 / 0.217 | ||
(tab_gabnfs <- tab_model(gab.nfs, show.p = F, show.se = TRUE))| ab land NFS | |||
|---|---|---|---|
| Predictors | Incidence Rate Ratios | std. Error | CI |
| (Intercept) | 374.47 | 125.81 | 204.76 – 711.73 |
| fc 5km | 0.95 | 0.01 | 0.93 – 0.97 |
| Observations | 12 | ||
| R2 Nagelkerke | 0.817 | ||
(tab_brc2nfs <- tab_model(rc2.nfs, show.p = F, show.se = TRUE))| brc abund nfs NA | |||
|---|---|---|---|
| Predictors | Estimates | std. Error | CI |
| (Intercept) | 0.82 | 0.34 | 0.07 – 1.58 |
| fc 3km | -0.02 | 0.01 | -0.05 – 0.00 |
| Observations | 12 | ||
| R2 / R2 adjusted | 0.267 / 0.194 | ||
FIGURES
Gamma hab
### Gamma ####
# r2(g.fs) # 0.55
## prediction g2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(g.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
# r2(g.nfs) # 0.73
## prediction g1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds2 <- predict(g.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds2$fit
new2$ic.up <- preds2$fit + preds2$se.fit*1.96
new2$ic.low <- preds2$fit - preds2$se.fit*1.96
## Plot gamma richness ~ fc %
(g.hab <- ggplot(data= new, aes(x= fc_5km, y = pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data= data_div, aes(x=fc_5km, y=gamma_FS), col= "#0072B2") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
geom_point(data= data_div, aes(x=fc_3km, y=gamma_NFS), col= "#D55E00") +
xlab(" ") + ylab("Gamma") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=20.5, x=38,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.55"))) +
annotate("text", y=19, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.73"))) + labs(tag="a"))Alpha hab
# r2(a.fs) # 0.40
## prediction a2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(a.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
# r2(a.nfs) # 0.63
## prediction a1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(a.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot Alpha richness ~ fc_3km
(a.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=alpha_FS), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=alpha_NFS), col= "#D55E00") +
xlab("") + ylab("Alpha") +
scale_color_manual(breaks=c("FS", "NFS"),
values=c('FS'='#0072B2', 'NFS'='#D55E00')) +
theme(axis.text=element_text(size=15),
axis.title=element_text(size=14,face="bold"),
legend.position = 'bottom') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=13, x=38,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.40")))+
annotate("text", y=11, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.63"))) +
labs(tag="b"))Beta hab
### Beta add ####
# r2(b.fs) # 0.58
## prediction b1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(b.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
# r2(b.nfs) # 0.65
## prediction b1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(b.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot Beta Diversity ~ fc %
(b.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new, aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=betad_FS), col= "#0072B2") +
geom_line(data= new2, aes(x= fc_3km, y= pred),
col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=betad_NFS), col= "#D55E00") +
xlab(" ") + ylab("Beta") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(0,max(data_div[2:7])+0.5) +
# scale_x_reverse() +
annotate("text", y=18, x=40,colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.58")))+
annotate("text", y=16, x=38, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.65"))) +
labs(tag="c"))g.hab + a.hab + b.habggsave("figures/plot.hab.jpeg",units="cm", width=20, height=8, dpi=300)Total abundance FS & NFS
## FS
# r2(gab.fs) # 0.44
## prediction gab1
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(gab.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
# r2(gab.nfs) # 0.82
## prediction gab2
new2 <- new
preds <- predict(gab.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot gamma richness ~ forest cover 5km
(gab.hab <- ggplot(new, aes(x=fc_5km, y=pred)) +
geom_line(linetype= "dashed", col= "#0072B2", aes(alpha= 0.1)) +
# geom_ribbon(aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=ab.land_FS), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_5km, y=pred), col= "#D55E00") +
geom_ribbon(data= new2, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.2, fill= "#D55E00") +
geom_point(data= data_div, aes(x=fc_5km, y=ab.land_NFS), col= "#D55E00") +
scale_x_continuous(limits= c(10,50)) +
xlab(" ") + ylab("Total abundance") +
theme(axis.text=element_text(size=15),axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
# scale_x_reverse() +
annotate("text", y=300, x=38, alpha= .5, colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.44"))) +
annotate("text", y=270, x=40, colour= "#D55E00",
label=expression(paste(~ R^2 ~ "= 0.82"))) +
labs(tag="d"))
### Beta RC FS & NFS
## FS
# ## NULL ## r2(rc1.fs)
## prediction rc2.2 and rc2
new <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=30))
preds <- predict(rc1.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
# r2(rc1.nfs) # 0.25
## prediction rc2
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc1.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## Plot RCbeta abundance ~ forest cover at 3
(rc.hab <- ggplot(data= new, aes(x=fc_3km, y=pred)) +
geom_line(data= new,aes(x=fc_3km, y=pred),
col= "#0072B2", linetype= "dashed", alpha= .2) +
# geom_ribbon(data= new, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.fsNA), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), col= "#D33E00", linetype= "dashed", alpha= .2) +
# geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D33E00") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.nfsNA), col= "#D33E00") +
xlab("Percent forest cover") + ylab(c(expression(bold("β"["RC"])))) +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),
legend.position = 'none') + ylim(-1,1) +
# scale_x_reverse() +
annotate("text", y=0.9, x=43, alpha= .5, colour= "#D33E00",
label=expression(paste(~ R^2 ~ "= 0.25"))) +
# annotate("text", y=0.7, x=43, alpha= .5, colour= "#0072B2",
# label=expression(paste(~ R^2 ~ "= 0.30"))) +
# annotate("text", y=0.9, x=43, colour= "#D33E00",
# label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.37"))) +
# annotate("text", y=0.8, x=43, colour= "#0072B2",
# label=expression(paste(beta, "(SE) = -0.06(0.03); " ~ R^2 ~ "= 0.33"))) +
labs(tag="e"))
Beta RC abund FS & NFS
# r2(rc2.fs) # 0.41
## prediction rc2
new <- data.frame(fc_5km = seq(min(data_div$fc_5km),max(data_div$fc_5km), length.out=50))
preds <- predict(rc2.fs, type="response", newdata=new, se.fit=T)
new$pred <- preds$fit
new$ic.up <- preds$fit + preds$se.fit*1.96
new$ic.low <- preds$fit - preds$se.fit*1.96
## NFS
# r2(rc2.nfs) # 0.21
## prediction rc1
new2 <- data.frame(fc_3km = seq(min(data_div$fc_3km),max(data_div$fc_3km), length.out=50))
preds <- predict(rc2.nfs, type="response", newdata=new2, se.fit=T)
new2$pred <- preds$fit
new2$ic.up <- preds$fit + preds$se.fit*1.96
new2$ic.low <- preds$fit - preds$se.fit*1.96
## RCbeta abund ~ FC%
(rc.ab.hab <- ggplot(data= new, aes(x=fc_5km, y=pred)) +
geom_line(data= new,aes(x=fc_5km, y=pred), col= "#0072B2") +
geom_ribbon(data= new, aes(x=fc_5km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#0072B2") +
geom_point(data=data_div, aes(x=fc_5km, y=brc.abund.fsNA), col= "#0072B2") +
geom_line(data= new2, aes(x=fc_3km, y=pred), linetype= "dashed", col= "#D55E00", alpha=.2) +
# geom_ribbon(data= new2, aes(x=fc_3km , ymax=ic.up, ymin=ic.low), alpha=0.1, fill= "#D55E00") +
geom_point(data=data_div, aes(x=fc_3km, y=brc.abund.nfsNA), col= "#D55E00") +
xlab(" ") + ylab(c(expression(bold("β"["RC-abundance"])))) +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14,face="bold"),legend.position = 'none') +
ylim(-1,1.1) +
# scale_x_reverse() +
annotate("text", y=0.9, x=45, colour= "#0072B2",
label=expression(paste(~ R^2 ~ "= 0.41"))) +
annotate("text", y=0.75, x=43, alpha= .5, colour= "#D55E00",
label=expression(paste( ~ R^2 ~ "= 0.21"))) +
labs(tag="f"))Final plot
(plot.hab <- (g.hab + a.hab + b.hab) /
(gab.hab + rc.hab + rc.ab.hab))ggsave("figures/figure2.jpeg",units="cm", width=28, height=15, dpi=300)